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Abstract 

This paper presents a new resolution strategy for multi-scale streamer discharge sim- 
ulations based on a second order time adaptive integration and space adaptive mul- 
tiresolution. A classical fluid model is used to describe plasma discharges, considering 
drift-diffusion equations and the computation of electric field. The proposed numerical 
method provides a time-space accuracy control of the solution, and thus, an effective 
accurate resolution independent of the fastest physical time scale. An important im- 
provement of the computational efficiency is achieved whenever the required time steps 
go beyond standard stability constraints associated with mesh size or source time scales 
for the resolution of the drift-diffusion equations, whereas the stability constraint related 
to the dielectric relaxation time scale is respected but with a second order precision. 
Numerical illustrations show that the strategy can be efficiently applied to simulate the 
propagation of highly nonlinear ionizing waves as streamer discharges, as well as highly 
multi-scale nanosecond repetitively pulsed discharges, describing consistently a broad 
spectrum of space and time scales as well as different physical scenarios for consecutive 
discharge/post-discharge phases, out of reach of standard non-adaptive methods. 
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1. Introduction 



In recent years, plasma discharges at atmospheric pressure have been studied for an 
increasing list of applications such as chemical and biological decontamination, aero- 
dynamic flow control and combustion [l|, In all these physical configurations, the 
discharges take usually the form of thin plasma filaments driven by highly nonlinear ion- 
izing waves, also called streamers. These ionizing waves occur as a consequence of the 
high electric field induced by the fast variations of the net charge density ahead of an 
electron avalanche with large amplification. The streamer discharge dynamics are mainly 
governed by the Courant, the effective ionization and the dielectric relaxation times scales 
[3(, which are usually of the order of I0 -14 — I0 _12 s, whereas the typical time scale of 
the discharge propagation in centimeter gaps, is about a few tens of nanoseconds. On 
the other hand, a large variation of space scales needs also to be taken into account, 
since the Debye length at atmospheric pressure can be as small as a few micrometers, 
while the inter-electrode gaps, where discharges propagate, are usually of the order of a 
few centimeters. As a result, the detailed physics of the discharges reveals an important 
time-space multi-scale character 0, Q ■ 

More complex applications include plasma assisted combustion or flow control, for 
which the enhancement of the gas flow chemistry or momentum transfer during typical 
time scales of the flow of 10~ 4 — I0 _3 s, is due to consecutive discharges generated by 
high frequency (in the kHz range) sinusoidal or pulsed applied voltages (a, Li- Therefore, 
during the post-discharge phases of the order of tens of microseconds, not only the time 
scales arc very different from those of the discharge phases of a few tens of nanoseconds, 
but a completely different physics is taking place. Then, to the rapid multi-scale con- 
figuration during discharges, we have to add other rather slower multi-scale phenomena 
in the post-discharge, such as recombination of charged species, heavy-species chemistry, 
diffusion, gas heating and convection. Therefore, it is very challenging to accurately 
simulate the physics of plasma/flow interactions due to the synergy effects between the 
consecutive discharge/post-discharge phases. 

In most numerical models of streamer discharges, the motion of electrons and ions is 
governed by drift-diffusion equations coupled with Poisson's equation. Early simulation 
studies were limited to simplified situations where the streamer is considered as a cylinder 
of constant radius [1, 0, in which the charged particle densities are assumed 
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to be constant along the radial extension of the streamer: the 1.5D model approach. 
In this model, the spatio-temporal evolution of the charged particle densities is solved 
only along one spatial dimension in the direction of propagation, whereas the electric 
field is calculated in two dimensions using the so-called disc method, based on a direct 
integration of analytical results. A 2D model for the electric field is indeed essential to 
properly calculate the electric field enhancement by the space charge in the streamer 
head. After the first 2D streamer simulations using the Poisson's equation resolution 
were performed Hp , m any studies have been carried out in 2D 

IE EE El EE EE EE El 



and 3D [H, HE 

Being aware of the complexity of fully coupled resolutions of these modeling equa- 
tions, a decoupling strategy is usually adopted, which considers an independent and 
successive numerical resolution of Poisson's equation with a fixed charge distribution, 
and of the drift-diffusion equations with a fixed electric field during each decoupling time 
step. These computations might be performed explicitly in time with standard first or 
even second order schemes 23J, |24j. In these cases, the time steps are usually limited 
for the sake of stability by the various characteristic times scales (Courant, ionization, 
dielectric relaxation), whereas the accuracy of simulations is assumed to be given by the 
resolution of the fastest physical time scale. In order to somehow overcome the dielectric 
relaxation limitation, some semi-implicit approaches were developed [IE El US' based 
on a predictive approximation of the space charge ahead in time during the electric field 
computation, even though the other time scale constraints remain. This gain of stability 
allows important improvements in terms of computational efficiency but the accuracy of 
simulations becomes rather difficult to quantify. 

Another performing technique to improve the efficiency of simulations considers an 
asynchronous explicit time integration of the drift-diffusion equations with self-adaptive 
local time-stepping, for which the local time steps are based either on local dynamic 
increments of the solution 
are the subject of several studies 



291 or on local Courant conditions [30( . These techniques 
33j] and are mainly conceived to avoid expensive 



31, 32 



computations whenever the whole system is unnecessarily advanced in time with a global 
time step prescribed by the fastest scale. Even though these methods yield efficient 
strategies, specially in terms of CPU time savings, with stable and flux-conserving time 
integrations, it is rather difficult to conduct an accuracy control on the resolution of 
the time dependent equations or on their coupling with the electric field resolution for 
plasma models. 

In this work, a numerical study is conducted in order to build a second order explicit 
in time decoupling scheme for the resolution of the electric field and the electron and 
ion densities. A lower order and embedded method is taken into account to dynamically 
compute the decoupling time steps that guarantee an accurate description with error 
control of the global physical coupling. At this stage, the only limiting time scale is the 
dielectric relaxation characteristic time for stability reasons. In a second level, the drift- 
diffusion equations are solved using a Strang second orde r op erator splitting scheme 



in order to guarantee the global order of the strategy [34|, |35[. This time integration 



scheme considers high order dedicated methods during each splitting time step, which is 
dynamically adapted by an error control procedure |36j |. In this way, even though there 
is a global advance in time given by the splitting time step, the latter is determined by 
the desired accuracy of the global physics, which is not necessarily related to the stability 
constraints associated with the mesh size or the fastest source time scales as demonstrated 
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in [37[ . As a consequence, this technique provides an error control procedure and stands 
as an alternative way to local stepping schemes to overcome time step limitations related 
to the reaction, diffusion and convection phenomena. 

Both the electric field and density resolutions are performed on an adapted mesh 
obtained by a spatial multiresolution method, based on Harten's pioneering work f38j ] 
and further developed in [39| , taking into account the spatial multi-scale features of these 
phenomena with steep spatial gradients. In particular, some grid adaptation techniques 
for 2D structured meshes were already used 23, 3,[20[ and extensions to 3D have been also 
proposed (20l . ITjj for streamer simulations. However, one of the main advantages of the 
multiresolution approach is that it is based on a wavelet representation technique and an 
error of the spatial approximation can be then mathematically estimated. Consequently, 
an effective error control is achieved for both the time and space resolution of the multi- 
scale phenomena under study. 

The performance of the method is first evaluated for a propagating streamer problem 
with the multi-scale features previously discussed, for which the various simulation pa- 
rameters are studied. Once the physical configuration is settled, a 1.5D streamer model 
is adopted in order to obtain an electric field resolution strategy based on direct com- 
putations and derived from analytical expressions, suitable for adapted finite volume 
discretizations 0. In a second step, a more complex physical configuration is consid- 
ered for the simulation of repetitively pulsed discharges, for which a time-space adaptive 
method is required to efficiently overcome some highly multi-scale features in order to 
fully describe the various physical phenomena. In this work, only a 1.5D model is consid- 
ered but extensions to higher dimensions is straightforward for instance with a Poisson's 
equation solver for adapted grids as it has been implemented in (23l. 13. |20| . However, 
in this paper we focus on the development and validation of new numerical methods 
for the resolution of the drift-diffusion equations and its coupling with the electric field 
computation, which arc independent of the dimension of the problem. Numerical illus- 
trations of multidimensional problems with the same time-space adaptive strategy with 
error control will be the subject of future work. 

The paper is organized as follows: in Section^ we present the physical configuration 
and the modeling equations. The numerical strategy is presented in Section [31 in which 
the second order adaptive time integration technique is detailed along with the resolution 
of drift-diffusion equations and the electric field, as well as the spatial multiresolution 
adaptive procedure. Numerical illustrations are summarized in Section [4] for two con- 
figurations given by single propagating and multi-pulsed discharges. We end with some 
concluding remarks and prospects on future developments and applications. 



2. Model formulation 

In this work, we consider positive streamer discharges in air at atmospheric pressure 
in a point-to-plane geometry, as shown in Figure [TJ The tip of the anode is placed 1 cm 
from the planar cathode and the radius of curvature of the anode is 324 /zm. The most 
common and effective model to study streamer dynamics is based on the following drift- 
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diffusion equations for electrons and ions, coupled with Poisson's equation (41. 42) : 



d t n e - <9 X 
d t n p 

d t n n - <9 X • n n v n - <9 X ■ (D n <9 x n n ) 



d x ■ ripVp - <9 X ■ (D p d x n p ) 



n a\v c \ - n ?7|v c | + n c n p /? cp + n„7, 
n c a\v c \ - n c n p f3 C p + n n n p f3 np , )■ (1) 

n e 7?|v | - n n n p f3 n p - n n j, 



E Q d^V = -q c (rip 



n e ), 



(2) 



where x S K d , m is the density of species i (e: electrons, p: positive ions, n: negative 
ions), V is the electric potential, Vj = /i;E (E being the electric field) is the drift velocity. 
Di and /ij, are the diffusion coefficient and the absolute value of mobility of the charged 
species i, q Q is the absolute value of an electron charge, and Sq is the permittivity of free 
space, a is the impact ionization coefficient, r/ stands for the electron attachment on 
neutral molecules, /3 ep and /3 np account respectively for the electron-positive ion and the 
negative-positive ion recombination, and 7 is the detachment coefficient. 
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Figure 1: Computational domain for the studied point-to-planc geometry. 



The electric field E and the potential V arc related by 

E = -d x V, (3) 

and thus, the Poisson's equation (|2|) becomes: 

e q <9 x ■ E = q e (n p - n n - n c ). (4) 

All the coefficients of the model are assumed to be functions of the local reduced 
electric field E/N gas , where E is the electric field magnitude and is the air neutral 
density. For test studies presented in this paper, the transport parameters for air are 



taken from [43|; detachment and attachment coefficients, respectively from |44j and )45l |: 
and other reaction rates, also from [4^|. Diffusion coefficients for ions are derived from 
mobilities using classical Einstein relations. 

In simulations of positive streamer discharges in air at atmospheric pressure without 
any prcionization, the photoionization term is crucial to produce seed charges in front of 
the streamer head and then to ensure the streamer propagation [zjj . However, in repeti- 
tive discharges, [46[ and recently (47l | have shown that even at low frequency, a significant 
amount of seed charges from previous discharges may be present in the inter-electrode 
gap. In this work, we have neglected the photoionization source term and considered 
discharge conditions with a prcionization background to ensure a stable propagation of 



the discharge without impacting the main discharge characteristics [15|, |46|, |18 
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3. Construction of the numerical strategy 



In this section, we introduce a new numerical technique for multi-scale streamer 
discharge simulations, based on a second order decoupled resolution of the electric field 
and the drift-diffusion equations for electrons and ions, with self-adaptive decoupling time 
steps with error control. The drift-diffusion equations are then solved using a dedicated 
Strang time operator splitting scheme for multi-scale phenomena. On the other hand, 
the electric field is computed based on a parallel computing method, specially conceived 
for the configuration under study in 1.5D geometry. Both resolutions are conducted on 
a dynamic adaptive mesh using spatial multircsolution transformation with error control 
of the spatial adapted representation. 

3.1. Second order adaptive time integration strategy 

Let us write the semi-discretized equations (HJ and ^ in the following way just for 
analysis purposes: 

d t ip = ${ip,<f>), \ , , 

o = $(,M), / W 

for t > t 0l where ip : R -)• R Nxm and : M -> R Nxd stand respectively for the 
spatial discretization of (n c ,n p ,n n ), i.e. m = 3, and of E over N points. Supposing 
that all functions are sufficiently diffcrcntiablc in all their variables and using the Taylor 
expansion of the true solution, one can write after some time At from initial time to, 

At 2 

i>(t + At) = Vo + At*(^ , <h) + — [fy* * + fy* M] t=t0 + 0{At 3 ), (6) 

with ip = 4'(t Q ), <t>o = 4>{to). 

A second order in time resolution of system ([5]) must then verify Q locally for each 
At. However, as it was stated before, solving simultaneously ([T]) and ^ (or or 
equivalently ([5]), involves important numerical difficulties, considering for instance the 
different nature of equations (JTJ) and ([2]) (or ([!])). Therefore, a decoupled approach is 
often used in which one aims at solving the drift-diffusion equations and the electric field 
independently. This amounts to solve 

dt^ = #(iM*). te}t Q ,t Q + At], (7) 

with fixed 4>* = </>(**), t* € [t ,t + At] and 4>(t ) = ip . 

The most common technique considers t* = to, that is, to previously compute the 
electric field at to from $(?/>o, </>o) = 0, and then solve (JT]) with (j>* = 4> . This can be 
interpreted as a standard first order operator splitting method that yields an approxima- 
tion of order 1, ipi(t), of the exact solution, ip(t), based on classical numerical analysis 
results obtained by confronting §Q with 

At 2 

Mto + At) = Vo + Am^o,^) + — [d^*} t=to + 0(At 3 ). (8) 

The same result follows for </>i(to + At) computed out of $>(ipi(to + At), 4>i(to + At)) = 
or equivalently, out of its explicit representation 0i(to + At) = T(^i(fo + At)), assuming 
a Lipschitz condition: 

||TW-T(^)||<i]|V-^ll- (9) 
6 



Considering now any t* G [t ,t +At] into ([7]), the only second order solution, (^(t), 02(t)), 
will be given by resolution of (O with 0* = 0i for t* = to + At/2, for which 



where 



At 2 

V> 2 (*o + A*) = + Ai*(Vo, 0i ) + [fy* 4-] 

2 z 



At 



0(At 3 ), 



(10) 



#(V>O,0i) = * Uo,0 to^- 



£>(At 2 ) 



*(^o,0o) + ^ 



Jt=t 



0(At 2 ), 



and hence, 



At 2 



Mto + At) = V^o + At*(<0o, 0o) + — [fy* * + d*0] t=to 



0(At 3 ); 



(11) 
(12) 

2 (t o + At) = T(V> 2 (t + At)). (13) 

Nevertheless, this second order approximation, ^(t), is based on the previous knowl- 
edge of 0i = 0(to + At/2), and thus, of ip(to + At/2). In order to overcome this difficulty, 

one can solve (J7J with 0* = 0i(to + At/2) = T(i/»i(to + At/2)), that is, computing first 
■0i(to+At/2) with the first order method. In particular, this does not change the previous 
order estimates as it follows from 



and 



V>(io + At) -M*o + At) = 



At 2 

2 
At 2 



d^dt 



0(At 3 



t=t 



£>(At 3 ) 



C(At 3 ). 



Taking into account both methods 

T-At 

- 'l 



M*o + At) 
0i (to + At) 



00 
00 



V> 2 (to + At) 
2 (t o + At) 



7? 



At ( IpQ 



(14) 



(15) 



we perform computations with a second order scheme 7^* , which uses an embedded and 
lower order scheme 7^ A *^ 2 , as it was previously detailed. An adaptive time step strategy 
is then implemented in order to control the accuracy of computations by tuning the 
duration of the decoupled resolution. It is based on a local error estimate, dynamically 
computed at the end of each decoupling time step At, given by 

||7f*(Vo,0o)*- 7^0,00)1 



C(At 2 ). 



Therefore, for a given accuracy tolerance 777-, 



\Tf\ik,<h) t -Tf\ih,ih) t \\<rrr 
7 



(16) 
(17) 



must be verified in order to accept the current computation with At, while the new time 
step is calculated by 



Ai ncw - At J r^AtZT 7^ _ rr Ml , , vll • (18) 



vr 

|7^*(Vo,0o)* -T^o^oY 

Several dedicated solvers can be then implemented for each subproblem (|TJ) and ^ 
while the theoretical error estimates of the decoupling schemes analyzed in this section 
remain valid. In this way, the independent choice of appropriate numerical schemes 
allows to strongly reduce the computational complexity of the global numerical strategy, 
and an error control procedure such as the one proposed in this work allows to effectively 
calibrate this decoupling within a prescribed accuracy tolerance. 

3.2. Resolution of the drift- diffusion equations 

We consider now the numerical resolution of the drift-diffusion equations ([T]), that 
can be written in the general form of a convection-reaction-diffusion system of equations: 



a t u-a x (F(u)+D(u)9 x u) = f(u), t>t , 
u(*o,x) = u (x), t = t , 



(19) 



where F, f : R ln — > R m and u : R x R d — > R m , with a tensor of order d x d x m as 
diffusion matrix D(u). In particular, u = (n e , n p , n n ) with m = 3 in this study. 

The system (fT9"]l corresponds to problem ([7]) for a fixed electric field, and it is solved 
during each decoupling time step At into T2 (or 71 ) scheme, using a Strang time operator 
scheme with dedicated high order time integrators on a dynamic adaptive mesh, based 
on a strategy introduced in [34|]. This resolution is briefly detailed in following sections. 

3.2.1. Time operator splitting 

An operator splitting procedure allows to consider dedicated solvers for the reaction 
part which is decoupled from other physical phenomena like convection, diffusion or both, 
for which there also exist dedicated numerical methods. These dedicated methods chosen 
for each subsystem are then responsible for dealing with the fastest scales associated with 
each one of them, in a separate manner, while the reconstruction of the global solution 
by the splitting scheme should guarantee an accurate description with error control of 
the global physical coupling, without being related to the stability constraints of the 
numerical resolution of each subsystem. 

Considering problem (|19[) and in order to remain consistent with the second order T2 
scheme, a second order Strang scheme is implemented [3] 

«S At °(u ) = n Au ' 2 V At -' 2 C A ^V At ^ 2 n Au ' 2 {n ), (20) 

where operators 1Z, T>, C indicate respectively the independent resolution of the reaction, 
diffusion and convection problems with a splitting time step, At s , taken inside the overall 
decoupling time step, At s < At. Usually, for propagating reaction waves where for 
instance, the speed of propagation is much slower than some of the chemical scales, the 
fastest scales are not directly related to the global physics of the phenomenon, and thus, 



larger splitting time steps might be considered [34l. |35| . Nevertheless, order reductions 



may then appear due to short-life transients associated with fast variables and in these 



cases, it has been proven in [50( that better performances are expected while ending the 
splitting scheme by operator 1Z or in a more general case, the part involving the fastest 
time scales of the phenomenon. 

The resolution of (|TT)|) should be precise enough to guarantee theoretical estimates 
given in Section [3. II Therefore, an adaptive splitting time step strategy, based on a local 
error estimate at the end of each splitting time step At s , is also implemented in order to 
control the accuracy of computations [37|]- In this context, a second, embedded and lower 
order Strang splitting method S Ats was developed by [36|, which allows to dynamically 
calculate a local error estimate that should verify 

||<S Ats (u ) -«S At =(u )|| « 0(At s 2 ) < r Mit , (21) 

in order to accept the current computation with At s , and thus, the new splitting time 
step is given by 




a '" " 1 ^ ^HuJ^Mr <° + At - 1 1 • (22) 

with 77 sp nt < ?yr and t = J^i At Si while t € ]t , t + At]. 

The choice of suitable time integration methods to numerically approximate TZ, T> and 
C during each At s is mandatory not only to guarantee the theoretical framework of the 
numerical analysis but also to take advantage of the particular features of each indepen- 
dent subproblem. A new operator splitting for reaction-diffusion systems was recently 



introduced [34j, [35j, which considers a high fifth order, ^4-stable, L-stable method like 
Radau5 [5l|, based on implicit Rungc-Kutta schemes for stiff ODEs, that solves with a 
local cell by cell approach the reaction term: a system of stiff ODEs without spatial cou- 
pling in a splitting context. For the diffusion problem, another high fourth order method 
like ROCK4 [HH is considered, which is based on explicit stabilized Runge-Kutta schemes 
that feature extended stability domains along the negative real axis. The ROCK4 solver 
is then very appropriate for diffusion problems because of the usual predominance of neg- 
ative real eigenvalues. Both methods incorporate adaptive time integration tools, similar 
to (|T5|) and (|2"2"j) . in order to control the accuracy of the integrations for given accuracy 
tolerances r?R ad au5 and ?7rock4, chosen such that ??Radau5 < Tfeplit and t7 RO ck4 < i? sp ut - In 
particular, in the case of multi-scale propagating waves, it can be proven that the local 
treatment plus the adaptive time stepping of the reaction solver allow to discriminate the 
cells of high reactive activity only present in the neighborhood of the localized wavefront, 



saving as a consequence a large quantity of integration time 35 1. 

An explicit high order in time and in space one step monotonicity preserving scheme 
OSMP [53J is used as convective scheme. It combines monotonicity preserving constraints 
for non-monotone data to avoid extrema clipping, with TVD features to prevent spurious 
oscillations around discontinuities or sharp spatial gradients. Classical CFL stability 
restrictions are though imposed inside each splitting time step At s for operator C A * S . 
The overall combination of an explicit treatment of the spatial phenomena as convection 
and diffusion, with a local implicit integration of stiff reaction implies important savings 
in computing time and memory resources [34{ . as well as an important reduction of 
computational complexity with respect to a fully implicit coupled resolution of problem 
(|T9"|) . On the other hand an explicit coupled treatment of (|T9"|) will have a very limited 
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efficiency for stiff problems unless more sophisticated strategies as the asynchronous local 
time-stepping techniques [28, 2i| 3(| are considered even though these schemes do not 



provide a precise measurement of the accuracy of the integration. 

Finally, the numerical errors of the splitting scheme arc effectively handled by an error 
control procedure which furthermore allows to determine the coupling time scales of the 
global phenomenon that can be several orders of magnitude slower than the fastest time 
scales of each subproblcm treated by each dedicated solver. In this way a decoupling of the 
time scale spectrum of the problem is achieved that leads to more efficient performances 
within a prescribed accuracy tolerance whenever this decomposition of scales is possible. 

3.2.2. Mesh refinement technique 

Regarding problem (|19p . we are concerned with the propagation of reacting wave- 
fronts, hence important reactive activity as well as steep spatial gradients are localized 
phenomena. This implies that if we consider the resolution of the reaction problem, a 
considerable amount of computing time is spent on nodes that are practically at (partial) 
equilibrium. Moreover, there is no need to represent these quasi-stationary regions with 
the same spatial discretization needed to describe the reaction front, so that the convec- 
tion and the diffusion problems might also be solved over a smaller number of nodes. 
An adapted mesh obtained by a multircsolution process which discriminates the various 
space scales of thephenomenon, turns out to be a very convenient solution to overcome 
these difficulties [13, Furthermore, in plasma applications, the resolution of Pois- 
son's equation takes usually ^80% of the computing time. Thus, important savings are 
achieved with a mesh adaptive technique, as a consequence of the strong reduction of 
cells. 

In practice, if one considers a set of nested spatial grids from the coarsest to the finest 
one, a multiresolution transformation allows to represent a discretized function as values 
on the coarsest grid plus a series of local estimates at all other levels of such nested grids 
[39| . These estimates correspond to the wavelet coefficients of a wavelet decomposition 
obtained by inter-level transformations, and retain the information on local regularity 
when going from a coarse to a finer grid. Hence, the main idea is to use the decay 
of the wavelet coefficients to obtain information on the local regularity of the solution: 
lower wavelet coefficients arc associated with locally regular spatial configurations and 
vice-versa. The basis of this strategy is presented in the following. For further details on 
adaptive multiresolution techniques, we refer to the books of j55l | and [56^ . 

3.2.3. Basis of a multiresolution representation 

To simplify the presentation let us consider nested finite volume discretizations of 
(TT9)) with only one component, m = 1. For I = 0, 1, • • • , L from the coarsest to the finest 
grid, we have then regular disjoint partitions (cells) (r2 7 ) 7S s ( of an open subset ft C M d , 
such that each f2~, 7 € Si, is the union of a finite number of cells tt^. (i G S/+i, and 
thus, Si and Si + i are consecutive embedded grids. We denote U; := (u 7 ) 7e g ; as the 
representation of u on the grid Si where m 7 represents the cell-average of u : M x M. d —> ffi 
in fi 7 , 

u 7 := l^l" 1 / u(i,x)dx. (23) 

The data at different levels of discretization are related by two inter-level transfor- 
mations which arc defined as follows: 
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1. The projection operator P/_ 1; which maps U; to U;_i. It is obtained through exact 
averages computed at the finer level by 

u 7 = |Q 7 |- 1 I mI*V- ( 24 ) 

|/i| = | 7 | + l, Q^cO, 

where I7I := I if 7 £ 5 1 ;. As far as grids are nested, this projection operator is eraci 
and unique [Hil ]. 

2. The prediction operator P{~ , which maps U;_i to an approximation U; of U;. 
There is an infinite number of choices to define -P/ -1 , but we impose at least two 
basic constraints [39^ : 

(a) The prediction is local, i.e., for a given f2^ depends on a set of values it 7 
in a finite stencil surrounding $7^, where = |7| + 1. 

(b) The prediction is consistent with the projection in the sense that 

|fi> 7 = V \%\u„; (25) 



1/^1— l7l+i,n M co 7 



With these operators, we define for each cell 51^ the prediction error or detail as the 
difference between the exact and predicted values: 

d^ := ?v - u M , (26) 

or in terms of inter-level operations: 

d^^-P^opM^. (27) 

We can then construct a detail vector D/ as shown in (39i | in order to get a one-to-one 
correspondence from expressions (|26p and (1251): 



U,<— ►(U,_i,D I ). (28) 

Hence, by iteration of this decomposition, we finally obtain a multi-scale representation 
of U L in terms of M L = (U , Di, D 2 , • ■ • , D L ): 

M : U L — > M L , (29) 

where the details computed with (|27|) stand for the wavelet coefficients in a wavelet basis. 

One of the main interests of carrying on such a wavelet decomposition is that this 
new representation defines a whole set of regularity estimators all over the spatial domain 
and thus, a data compression might be achieved by deleting cells whose detail verifies 

|<y<e ; , l = \ft\, e l = 2^ l -Vr]MR, (30) 

where t)mr is a threshold value defined for the finest level L (38j . 

An important theoretical result is that if we denote by V£ := (v")xeS L j tnc solution 
fully computed on the finest grid, and denote by UJ, the solution reconstructed on the 
finest grid that used adaptive multircsolution (keeping in mind that the time integration 
was really performed on a compressed representation of U n ); then, for a fixed time 
T = nAt, it can be shown [38ll39l| that the approximation error made by using this space 
adaptive representation is proportional to the threshold value ?7mr: 

\\XJZ-Vl\\ L 2 cxm ?M R. (31) 
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3.3. Computation of the electric field 

In this part, we are concerned with the resolution of the electric field according to 



the 1~2 (or 7i ) scheme at some fixed time for a given distribution of charges (n e 



nn), 



considering a 1.5D model. This computation is also performed on the adapted mesh 
obtained by the previous multircsolution analysis. 



3.3.1. Discretization of the computational domain 

According to Figure [TJ the computational domain is limited by a planar cathode 
at x — and the tip of a hyperbolic anode at x — L x . The anode is not included 
in the domain. We consider streamers of fixed radius R s along the axis of symmetry. 
The computational domain is divided into n x cells of different size corresponding to the 
multircsolution adapted mesh, with faces x\, where i 6 [0,n x ] and cell centers x J c , where 
j G [l,7ia;]. The face x® corresponds to the position of the cathode and x™* corresponds 
to the position of the tip of the anode. Therefore for each cell x l c , there is its left face 



Figure E}. 



and its right face x\. For each cell x 3 c we define a width Wj 



Cj (see 



Figure 2: Definition of the grid: the cell centers are located at x J c , 
The domain is bounded by faces (cathode) and x™* (tip of the anode). 



3.3.2. Resolution of the electric field in a 1.5D model 

To determine the electric field during the propagation of the streamer, the space 
charge of the streamer is considered as a set of finite cylinders of width Wj , bounded by 
cell faces x\~ and x\ . As the computational domain is bounded by conducting electrodes 
of fixed potential, each volume charge pj creates an infinite series of image charges [H, Q . 
Then the principle of superposition is used to sum individual contributions from all the 
cylindrical space charges in the domain, their image charges, and the Laplacian electric 
field (computed based on classical results 57|). An advantage of this approach dwells in 
the fact that the electric field contributions from individual cylinders can be expressed 
analytically in a simple form and the determination of the electric field in each point of 
the domain can be performed in parallel. 

In the configurations we have studied, the cathode is grounded whereas an electric 
voltage is applied on the anode. These boundary conditions are taken into account by 
the Laplacian electric field and by including a series of image charges of the charges in 
the gap. It is important to note that the computation of the Laplacian electric field takes 
into account the real geometry of electrodes as shown in Figured] However, in this work, 
to simplify the computation of image charges we have assumed that both electrodes are 
planar. For a volume charge pj centered at x J c , there exist image charges of the first 
order with charge — pj at x = 2L X — x J c mirrored through the anode, sec Figure [3^, and 
at x = —x 3 c mirrored through the cathode, sec Figure [3}d. And for each of these image 
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charges there exist higher order image charges of opposite signs and so forth. All the 
image charges of pj up to order three are depicted in Figure [3};. 




Figure 3: Image charges up to the third order: (a) charge pj is first mirrored behind the anode (x = L x ), 
(b) charge pj is first mirrored behind the cathode (x = 0), (c) charge pj and its images. 



Integrating the generalized Coulomb's law |58j and using the principle of superposi- 
tion, we find that the cylinder charges of cells j e [1, n x ] of width Wj, radius i? s , charged 
with densities pj (see Figure 0]), and the Laplacian electric field El{x\) at x\ [57[, create 
the electric field E at position x\ as follows: 



E{x\)=E L {x\) + Y,s 



Pj 1 



3 = 1 



2e n 



Wj + 2hi 



h 2 itj + R 2 s + yj(hi tj 



\ 



ffi 



(32) 



where 



''■,/ 



x\ ~ xj for i > j , 

j_i and s — 

x { — x\ for i < j 



-1 for i > j 
-1 for i < j 



hi i E M) 
i 




X 



Figure 4: Charged cylinder considered to compute the electric field in the 1.5D model. 



The positive sign of s accounts for the electric field calculated on the right from the 
position of the charged cylinder and vice-versa. The same formula applies for the image 
charges, but an appropriate sign of the charge has to be carefully taken into account 
according to Figure [3] In particular, in a shared memory computing environment, a 
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straightforward parallelization is accomplished for equation (|32j) , in which each core solves 
successively the electric field on one single position x\ , and where neither synchronization 
stages nor data exchange are needed among nodes. 

Note that for R s — > oo (infinite plane charges), equation (j3"2")l yields the exact electric 
field for a planar front: 

For finite radius i? s the solution (j3"2"|) is valid only on the axis of the discharge, but when 
applied to a discharge of a small radius, the electric field will vary only negligibly over 
the cross section of the discharge. This approach is expected to be more accurate for any 
finite radius than any discretization of Poisson's equation Q . 



4. Numerical results 



In this section, we present some numerical illustrations of the proposed numerical 
strategy for the simulations of positive streamers using a 1.5D model in a point-to- 
planc geometry. First, we will consider a discharge propagation with constant applied 
voltage for which different features of the numerical strategy are discussed, e.g., error 
estimates, data compression values and computing time, in order to properly choose 
the simulation parameters. Then, the potential of the method is fully exploited for a 
more complex configuration of repetitive discharges generated by high frequency pulsed 
applied voltages, followed by a long time scale relaxation, for which a complete physical 
description of the discharge and the post-discharge phases is achieved. 

4-1- Propagation of a positive streamer with constant applied voltage 

We consider a point-to-plane geometry with a 1 cm gap between the tip of the elec- 
trode and the plane, and a constant applied voltage of 13 kV at x = L x . For the following 
simulations, the discharge is initiated by placing a neutral plasma cloud with a Gaussian 
distribution close to the tip of the anode. The initial distributions of electrons and ions 
are then given by 

n e,p( x )\ t =o = "m M exp (-(a; - c) 2 /w 2 ) + n Q , n n (x)\ t=0 = 0, (33) 

where w = 0.027cm, c = 1cm, n max = 10 14 cm~ 3 , and with a preionization of no = 
10 8 cm" 3 . There are no negative ions as initial condition. The streamer radius is set to 
R s = 0.05 cm to have a typical electric field magnitude in the streamer head of 120 kV/cm 
[59( . Homogeneous Neumann boundary conditions were considered for the drift-diffusion 
equations. 

Two instances of the discharge propagation are shown in Figure [SJ for 12 nested 
grids equivalent to 4096 cells on the finest grid, L = 12, and for accuracy tolerances of 
77-7- = Tfepiit = t]mr = 10~ 4 ; the spatial refinement takes place only where it is required. 
Fine tolerances were chosen in all cases for the solvers, ?7Radau5 = ?7rocK4 = 10~ 7 , to 
guarantee accurate integrations. For all the simulation cases, the detail in each cell 
is taken as the maximum of the details computed according to (|27[) for each variable, 
where the prediction operator is a polynomial interpolation of degree 2, performed on 
normalized log u of the density variables in order to properly discriminate the streamer 
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heads from the highly ionized plasma channel; this logarithmic scale guarantees a correct 
spatial representation of the phenomenon as seen in Figure [5] for the density profiles. 

In order to perform an analysis of the numerical results, we define the reference 
solution as a fine resolution with the 1~2 scheme that considers a fixed decoupling time 
step, At = 10~ 14 s and a uniform grid of 4096 cells. For this reference solution, the 
memory requirements are acceptable and the simulation is still feasible, but it requires 
about 14 days of real simulation time on an AMD Opteron 6136 Processor cluster, while 
running the electric field computation in parallel on 16 CPU cores. In this case, the 
computation of the electric field, based on a direct integration of individual contributions 
of the charged cylinders, represents 80% of total CPU time per time step (about 3.2 s). 

First of all, we must verify the previous order estimates for the 71 and 1~2 schemes 
given in Section [3. II We consider as initial condition the reference solution at t = 10 ns. 
In order to only evaluate errors coming from the decoupling techniques, T\ and T2, we 
consider a fine splitting time step, At s = 10~ 14 s, to solve the drift-diffusion problem 
((T|) and a uniform grid; then, we solve ([5]) with both schemes for several decoupling 
time steps Ai, , and calculate the normalized L 2 error between the first /second order and 
reference solutions after time t = 2 1Q At s = 1.024 x 10~ n s. Figure [5] shows results with 
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Figure 6: Normalized L 2 errors between the reference and the 71 (first order) and 72 (second order) 
solutions for several decoupling time steps At on a uniform grid of 4096 cells. Top: electron (left) and 
positive ions (right); and bottom: negative ions. 



Atj = 2 l At s , where i <E [1, 10], which clearly verify first and second order in time for the 
7i and T2 schemes, respectively, and prove important gains in accuracy for same time 
steps. For instance, for At < 10 -12 s the second order scheme provides solutions with L 2 
errors at least 100 times lower than those obtained with the first order method. 

Figure [7] shows the time evolution of the normalized L 2 error for each variable between 
the time-space adapted and reference solutions for several tolerances, 777- = 7? S piit = 
VMR = 10~ 4 , 10 -3 , and 10 -2 . These are rather approximations of the error since the 
reference and adapted solutions are not evaluated exactly at the same time, and therefore, 
they are often slightly shifted of about ^10~ 14 — 10~ 13 s. In these tests, the decoupling 
time steps At were limited by the dielectric relaxation time step, AtoR, after noticing 
an important amount of rejections of computed time steps according to (|18p . whenever 
At >1.5 x AtnR- Otherwise, At is dynamically chosen in order to locally satisfy the 
required accuracy 777-, but it does not show important variations considering the self- 
similar propagating phenomenon. 

In Figure [51 we can see the corresponding adapted grid to each previous configuration 
with different tolerances. The representation of the electric field and the densities shows 
that for rj j- = ?7spiit = ?7mr = 10~ 2 , the streamer front propagates faster than in the 
reference case, with a slightly higher peak of the electric field in the front. On the other 
hand, for rjj- — r/ sp ut — tjur < 10 -3 , we observe a quite good agreement between the 
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Figure 7: Time evolution of the normalized L 2 errors between the reference and adapted solutions 
with r) = t]f = r7 S pii t = r^R = 10~ 4 , 10 — 3 , and 10~ 2 , and 4096 cells corresponding to the finest 
discretization. Top: electron (left) and positive ions (right); and bottom: negative ions. 



adapted and reference resolutions. 

We consider now an accurate enough resolution with 777- = ?7 8 plit = ')mr = 10 -4 and 
investigate the influence of the number of grids, that is, the finest spatial discretization at 
level L that should be taken into account. Figure |H] shows the adapted grids for L = 10, 
11 and 12, respectively equivalent to 1024, 2048 and 4096 cells in the finest grid; and a 
close-up of the corresponding electric fields in the discharge head at t = 8 ns. We see 
that for this level of tolerances, the streamer front propagates slightly slower than the 
reference case for L = 10, whereas L = 11 gives already good resolutions compared with 
the reference solution and with L = 12. In particular, higher values of L would need 
lower tolerances in order to retain regions at the finest level; this is already the case for 
L = 13 (equivalent to 8192 cells). Therefore, L = 11 with 2048 cells at the finest level 
seems to be an appropriate choice for this level of accuracy. 

Table Q] summarizes the number of cells in the adapted grid (# AG) at time t = 8ns, 
and the corresponding data compression (DC) defined as the percentage of active cells 
with respect to the equivalent number of cells for the finest discretization, in this case 
2048 for L = 11. For this propagating case, the data compression remains of the same 
order during the time simulation interval. The CPU computing times correspond to a 
time domain of study of t £ [0, 10] ns computed by one sole CPU core. If we consider 
for example total computing time for L = 11 and tolerances 777- = r) ap iit = f)MR = 10~ 4 , 
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Figure 8: Top: adapted grids (left) and electric fields (right) at t = 8 ns with 4096 cells corresponding 
to the finest discretization, and r\ = rjj- = r/ sp n t = tjmr = 10~ 4 , 10 -3 , and 10 — 2 . Bottom: zoom on the 
electron distributions (left) and the electric field (right) with the same parameters, and the reference 
solution. 




0.2 0.4 0.6 0.8 1 0.5 0.51 0.52 0.53 0.54 0.55 

x [cm] x [cm] 

Figure 9: Adapted grids (left) and electric fields (right) at t = 8 ns, for several finest spatial discretization 
L = 10, 11 and 12, ryj- = »; S pii t = t?mr = 10 — 4 , and the reference solution. 



it is ~44 times less expensive with respect to a resolution on a uniform grid with 2048 
cells and ijj- = %p\it = 10~ 4 (CPU time of 8552 s). This is quite reasonable, taking into 
account that the computing time for the electric field resolution is proportional to at 
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least 0(N 2 ) for N computing cells, after ([32]). 



Table 1: Number of cells in the adapted grid (#AG) and data compression {DC) at time t = 8ns, CPU 
computing time for t € [0, 10] ns, L = 11, and several tolerances r] = r)j- = /y sp iit = *?mr- 



V 


#AG 


DC% 


CPU(s) 


icr 6 


724 


35.35 


1360 


icr 5 


421 


20.56 


517 


10" 4 


263 


12.84 


193 


icr 3 


138 


6.74 


66 


icr 2 


70 


3.42 
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In conclusion, in this section we have shown that the numerical strategy developed can 
be efficiently applied to simulate the propagation of highly nonlinear ionizing waves as 
streamer discharges. An important reduction of computing time results from significant 
data compression with still accurate resolutions. In addition, this study allows to properly 
tune the various simulation parameters in order to guarantee a fine resolution of more 
complex configurations, based on the time-space accuracy control capabilities of the 
method. 

4-2. Simulation of multi-pulsed discharges 

In this section, we analyze the performance of the proposed numerical strategy on 
the simulation of nanosecond repetitively pulsed discharges 0, [6^1 ■ The applied voltage 
profile for this type of discharges is a high voltage pulse followed by a zero voltage 
relaxation phase. The typical pulse duration is ~10 -8 s, while the relaxation phase takes 
over ~10~ 4 s. The detailed experimental study of these discharges in air has shown that 
the cumulative effect of repeated pulsing achieves a steady-state behavior (gfj. In the 
following illustrations, we choose a pulse duration of T p = 15 ns, which is approximately 
equal to the time that is needed for the discharge to cross the inter-electrode gap. The 
rise time considers the time needed to go from zero to the maximum voltage and it is 
set to T r = 2 ns. The pulse repetition period is set to Tp = 10 -4 s, equal to 10 kHz of 
repetition frequency, a typical value used in experiments Q . We model the voltage pulse 
P by using sigmoid functions 

P(t, s,r,p) = l- a(-t, -s, r) - a(t, s+p 7 r), (34) 

with 

a(t,s,r) = - ; — — rrr-, (35) 

v ' ' ; l + exp(-8(i-s)/r)' v ' 

for time t, where s indicates when the pulse starts; r is the rise time; and p is the pulse 
duration; t, s, r, p G [0, Tp]. With a maximum applied voltage Vm ax , the applied voltage 
V{t) is computed by 



V(t) = V max • P It — 



T P 



Tp,T T , T r , T p . (36) 



In repetitively pulse discharges at atmospheric pressure and 300 K, as discussed in [46l . l47j . 
electrons attach rapidly to O2 molecules during the interpulse to form negative ions 
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(characteristic time scale of 20 ns). Then, the rate of the plasma decay is determined 
by ion- ion recombination [45, 46, 47 1. When the next voltage pulse is applied, electrons 



are detached with a rate taken from [44j. Therefore, as initial condition we assume a 
distribution similar to the end of the intcrpulse phase with a homogeneous prcionization 
consisting of positive and negative ions with a density of 10 9 cm -3 . For electrons, we 
consider a low homogeneous background of 10 1 cm -3 . This small amount of electrons as 
initial condition has a negligible influence on the results. 

We set the tolerances to rjj- = 7? S piit = i]mr = 10~ 4 and consider L = 11 grid levels, 
equivalent to 2048 cells in the finest grid. As in the previous configuration, homoge- 
neous Neumann boundary conditions were considered for the drift-diffusion equations. 
Figure \TU\ shows the time evolution of the decoupling time steps and the applied voltage 
for the first six pulses, even though simulation was performed for 100 pulses, that is 
t e [0, 10 -2 ] s. This simulation took over 8h44m while running the electric field compu- 
tation in parallel on 6 CPU cores of the same AMD Optcron 6136 Processor cluster; this 
gives an average of 5.24 minutes per pulse period. Figure [TU] shows also the fourth pulse 
for which the steady-state of the periodic phenomenon was already reached and almost 
the same numerical performance is reproduced during the rest of computations. The time 
steps are about ^10~ n s during pulses, then increase from ^10~ 12 s up to about ^10~ 6 s 
during a period ^6000 times longer, for which standard stability constraints are widely 
overcome according to the required accuracy tolerance. Solving this problem for such 
different scales with a constant time step is out of question and even a standard strategy 
that considers the minimum of all time scales would limit considerably the efficiency of 
the method as it is shown in the representation. In this particular case, the dielectric 
relaxation is the governing time scale during the discharge as in the previous case with 
constant applied voltage, whereas the post-discharge phase is alternatively ruled by dif- 
fusive or convective CFL, or by ionization time scale, with all security factors and CFL 
conditions set to one in Figure IT0T 
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Figure 10: Time evolution of the applied voltage and the decoupling time steps At for a multi-pulse 
simulation for the first 6 pulses (left) and for the 4th one (right) with its subsequent relaxation. Rejected 
time steps are marked with black crosses, while the minimum time scale corresponds to the blue line. 



The computation is initialized with a time step included in the pulse duration. Nev- 
ertheless, after each relaxation phase, since the new time step is computed based on the 
previous one according to Q18f>. this new time step will surely skip the next pulse. In 
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order to avoid this, each time we get into a new period, that is when \t/Tp\ changes, 
we initialize the time step with At = 0.5T r = 1 ns. This time step is obviously rejected 
as seen in Figure 1101 as well as the next ones, until we are able to retrieve the right 
dynamics of the phenomenon for the required accuracy tolerance. No other intervention 
is needed neither for modeling parameters nor for numerical solvers in order to automat- 
ically adapt the time step needed to describe the various time scales of the phenomenon 
within a prescribed accuracy. 

Figure [TT] represents the time evolution of the data compression which ranges from 
~2% up to ~16% during each pulse period. Regarding only the electric field resolution 
with the same time integration strategy, a grid adaptation technique involves resolutions 
^39 to ^2500 times faster, based on a really rough estimate for 0{N 2 ) operations. 
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Figure 11: Time evolution of the data compression for a multi-pulse simulation for the first six pulses 
(left) and for the fourth one (right) with its subsequent relaxation. 



Figure H2] presents the discharge dynamics for the first period. First, we observe at 
t = 10 ns after the beginning of the pulse, the propagation of a positive streamer in the 
gap. In Section 14.11 a preionization of positive ions and electrons was used to ensure 
the positive streamer propagation. In this section, seed electrons ahead of the streamer 
front are created as the front propagates by detachment of negative ions initially present. 
We note that at 15 ns, which corresponds to almost the end of the plateau before the 
decrease of the applied voltage, the discharge has crossed ~0.75 cm of the 1 cm gap. As 
a consequence, during the voltage decrease and at the beginning of the relaxation phase 
where the applied voltage is zero, there is a remaining space charge and steep gradients 
of charged species densities in the gap. Then for t = 50 ns, Figure [12] shows that the 
electric field in the discharge is almost equal to zero except in a small area where steep 
gradients of the electric field are observed but with peak values of only 30V/cm. Wc 
have checked that this area corresponds to the location of the streamer head at the end 
of its propagation as it is seen in the representation. We note that in the post-discharge, 
electrons are attaching and then at t = 50 ns, the density of positive ions is almost equal 
to the density of negative ions in the whole gap. At t = 99972 ns, the densities of charged 
species have significantly decreased due to charged species recombination. However, it is 
interesting to note that the location of the previous streamer head can still be observed 
at the same location as at t = 50 ns, but with much smaller gradients of charged species 
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densities and a very small electric field. This final state is the initial condition of the 
second pulse with a non-uniform axial preionization with positive and negative ions and 
a much smaller density of electrons. 

After a few repetitive pulses, we have observed that the discharge dynamics reached 
a steady-state behavior as observed in the experiments. To show the characteristics of 
the discharge when the steady-state is reached, Figure [T3] shows the discharge dynamics 
of the 100th period. The sequence of images is the same as in Figure [TSJ At the end of 
the 99th pulse, we have observed that the axial distribution of charged species in the gap 
is uniform and that the level of preionization is 5 x 10 10 cm -3 positive and negative ions 
and 10 4 cm -3 electrons. We note that 10 ns after the beginning of the 100th pulse the 
propagation of the discharge is faster than for the first pulse. This faster propagation is 
mostly due to the higher preionization level of positive and negative ions in the gap in 
comparison of the first voltage pulse. We observe that for the 100th pulse, 15 ns after the 
beginning of the pulse the discharge has almost completely crossed the inter-electrode 
gap and then during the relaxation phase, there is no remaining space charge in the whole 
gap. Consequently, 50 ns after the beginning of the 100th pulse, axial distributions of 
all charged species are uniform. As already observed for the first pulse, at 50 ns after 
the beginning of the voltage pulse most electrons have attached and then, the density 
of positive ions is almost equal to the density of negative ions in the whole gap. We see 
that the corresponding electric field distribution is not uniform at 50 ns, but no steep 
gradients are observed as for the first voltage pulse. At t = 9999998 ns, that is to say at 
the end of the 100th period, we note that a very low electric field is obtained in the gap. 
An axially uniform distribution of charges is obtained with 5 x 10 10 cm~ 3 for positive and 
negative ions and 10 4 cm -3 for electrons, which was the initial condition of the 100th 
pulse. This demonstrates the existence of a steady-state behavior of these nanosecond 
repetitively pulsed discharges. 
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Figure 12: First period of pulsed discharges. Top: propagation of the discharge in the domain at t = 10 ns 
after the beginning of the pulse (left); and at t = 15ns (right). Bottom: relaxation on the short time 
scale t = 50ns; and end of the relaxation phase after t = 99972 ns (right). 
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Figure 13: Steady-state of pulsed discharges (last period). Top: propagation of the discharge in the 
domain at t = 9900010 ns after the beginning of the pulse (left); and at t = 9900015 ns (right). Bottom: 
relaxation on short time scale t = 9900050 ns; and end of the relaxation phase t = 9999998 ns (right). 
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5. Conclusions 

The present work proposes a new numerical strategy for multi-scale streamer simu- 
lations. It is based on an adaptive second order time integration strategy that allows 
to discriminate time scales-related features of the phenomena, given a required level of 
accuracy of computations. Compared with a standard procedure for which accuracy is 
guaranteed by considering time steps of the order of the fastest scale, the control er- 
ror approach implies on the one hand, an effective accurate resolution independent of 
the fastest physical time scale, and on the other hand, an important improvement of 
computational efficiency whenever the required time steps go beyond standard stability 
constraints. The latter is a direct consequence of the self-adaptive time step strategy 
for the resolution of the drift-diffusion equations which considers splitting time steps not 
limited by stability constraints for reaction, diffusion and convection phenomena. So 
far, the global decoupling time steps arc limited by the dielectric relaxation stability 
constraint but with a second order accuracy. Nevertheless, we have also demonstrated 
that the decoupling time steps are rather chosen based on an accuracy criterion. Besides, 
if a technique such as a semi-implicit approach is implemented, the same ideas of the 
proposed adaptive strategy remain valid. 

An adaptive multiresolution technique was also proposed in order to provide error 
control of the spatial adapted representation. The numerical results have proven a natural 
coupling between time and space accuracy requirements and how the set of time-space 
accuracy tolerances tunes the precise description of the overall time-space multi-scale 
phenomenon. As a consequence, the numerical results for multi-pulsed discharge config- 
urations prove that this kind of multi-scale phenomena, previously out of reach, can be 
successfully simulated with conventional computing resources by this time-space adap- 
tive strategy. And they also show that a consistent physical description is achieved for a 
broad spectrum of space and time scales as well as different physical scenarios. 

In this work, we focused on a 1.5D model in order to evaluate the numerical perfor- 
mance of the strategy. However, the dimension of the problem will only have an influence 
on the computational efficiency measurements but not on any space-time accuracy or sta- 
bility aspects. At this stage of development, the same numerical strategy can be coupled 
with a multi-dimensional Poisson's equation solver, even for adapted grid configurations 
as developed recently in [2^, 0, H3|- Finally, an important amount of work is still in 
progress concerning programming features such as data structures, optimized routines 
and parallelization strategies. On the other hand, numerical analysis of theoretical as- 
pects is also underway to extend and further improve the proposed numerical strategy. 
These issues constitute particular topics of our current research. 
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